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We study the high density region of QCD within an effective model obtained in the frame of the hopping 
parameter expansion and choosing Polyakov type of loops as the main dynamical variables representing the 
fermionic matter. To get a first idea of the phase structure, the model is analyzed in strong coupling expansion 
and using a mean field approximation. In numerical simulations, the model still shows the so-called sign prob- 
lem, a difficulty peculiar to non-zero chemical potential, but it permits the development of algorithms which 
ensure a good overlap of the Monte Carlo ensemble with the true one. We review the main features of the model 
and present calculations concerning the dependence of various observables on the chemical potential and on the 
temperature, in particular of the charge density and the diquark susceptibility, which may be used to character- 
ize the various phases expected at high baryonic density. We obtain in this way information about the phase 
structure of the model and the corresponding phase transitions and cross over regions, which can be considered 
as hints for the behaviour of non-zero density QCD. 



PACS numbers: 11.15.Ha, 12.38.Gc, 12.38.Aw 

I. INTRODUCTION 

The exploration of the phase diagram of matter at non-zero 
baryon density is a challenging and interesting problem. In 
particular, it has been emphasized that quark matter at ex- 
tremely high density may behave as a color superconductor 
(see Ref.llll for a recent review on the subject and references 
therein). Moreover, it is also expected that the phase diagram 
in the temperature-density plane shows multiple phases sepa- 
rated by various critical lines and, except for the high T, small 
/i region, not much is known about their exact position and na- 
ture. 

Lattice gauge theory calculations in various implementa- 
tions that try to evade the sign problem generated by the non- 
zero chemical potential have been mostly performed at small 
baryon density and high temperature, where they agree rea- 
sonably well with each other Here there is good evidence 
for the presence of a crossover instead of a sharp deconfin- 
ing transition. At large ji (baryon density), however, there are 
only few numerical results which need to be corroborated by 
using different methods. See [jj] for a review. 

The aim of this work is to understand the phase structure of 
high density, strongly interacting matter. Most work on QCD 
at non-zero density proceeds from the jj. = 0, T ^ Tc region 
and attempts to go as far as possible in the /i > domain. As 
an alternative one may consider the possibility to start from 
the large ft, domain and try to reach the region of interest from 
above. In the spirit of the /i = quenched approximation a 
'non-zero density quenched approximation' for /i > based 
on the double limit 7\/^oo,yU^cxD,^ = exp (/i — In M) : 
fixed iH 01 has been considered. This implements a static, 
charged background, which influences the gluonic dynamics 
ii, 5]. The present model i6t| represents a systematic exten- 
sion of the above considerations: the gluonic vacuum is en- 
riched by the effects of dynamical quarks of large (but not 



infinite) mass, providing a large net baryonic charge. In |01 
and in the present paper we explore the phase structure of the 
model, as a first step in understanding the properties of such a 
background. 

This model can be derived as a 1/M expansion of QCD at 
large ^ around the unphysical limit of infinitely heavy quarks. 
However, it is more realistic to understand it as an approxima- 
tion whose justification relies on the predominant role of the 
gluonic dynamics. We want to understand how this dynamics 
is influenced by the presence of charged matter. This would 
allow, among other things, to study the effect of dense, heavier 
background baryonic charges on light quarks and hadrons. 

The main ingredient of the model are Polyakov-type loops, 
capturing the effect of heavy quarks with low mobility. The 
model still has a sign problem, but being based on the vari- 
ables which are especially sensitive to the physics of dense 
baryonic matter it allows for reweighting algorithms which 
ensure a good overlap of the Monte Carlo ensemble with the 
true one. 

The paper is organized as follows. In Sec|ll]we study the 
high density region of QCD within an effective model ob- 
tained by an expansion in the hopping parameter k of the 
fermionic determinant up to next-to leading order, k^. In 
Sec|lII]the model is analyzed using first a strong coupling ex- 
pansion and then a mean field approximation just to get a first 
idea of the phase diagram and to compare with numerical sim- 
ulations. 

SeclIVIshows results of the numerical simulations. Here the 
model shows the so-called sign problem but due to the factor- 
ization of the fermionic determinant it permits to develop very 
efficient local algorithms and achieve large statistics. The de- 
pendence of various observables on the chemical potential and 
the temperature is studied and we show a tentative phase dia- 
gram at large mass and high baryon density. Conclusions and 
outlook are given in Sec. |V] 
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II. QCD AT LARGE CHEMICAL POTENTIAL 



A. QCD at non-zero /x 

In this study we use the grand canonical formulation of 
QCD, i.e., we introduce the chemical potential /i as a (bare) 
parameter. The QCD grand canonical partition function with 
Wilson fermions at /i > is: 

Z(/3,K,7G,7F,/i) 

[Du] c 7^., {U}) , (2.1) 

5g(/3,7g,{C/}) 

/3 



Nr. 



ReTi 



7i., /i, {[/}) = Det W^K 7F, /i, {[/}) , 

3 

Wff = 6ff, [1 - Kf ^ (r+, a, + r_, t* u*) 

-KfjF{e^f r+4 U4T4 + C -^f r_4 jj;)] 
r±,, = 1 ± 7^, 7^ = 7*, 7^ = 1 , 

1 1 



(2.2) 
(2.3) 



2(A/ + 3 + 7F cosh/x) 2(Mo + 3 + 7f) ' 

where we have specialized Sg for Wilson's plaquette (P) ac- 
tion and used a certain definition of the Wilson term in W. 
Here AI is the 'bare mass', AIq the bare mass at /i = 0, / is 
the flavor index, denote the link variables and T"^ lattice 
translations. For the sake of generality and the discussion in 
section III.B we also introduced coupling anisotropics 7^, jp 
which however will be set to 1 elsewhere. All quantities are 
understood in units of the (spatial) lattice spacing a unless ex- 
plicitly specified otherwise. The exponential prescription for 
/i ensures canceling of divergences in the small a limit (iSfl. A 
non-zero physical temperature T is introduced as 



aT 



Iphys 
Nr 



(2.4) 



where jphys is the physical cutoff anisotropy defined by an 
appropriate renormalization of the coupling anisotropics i^, 
and Nr the 'length' of the (periodic) temporal lattice size. 
The fermionic coupling matrix W fulfills: 

75^^(^)75 = Wi-^i)*, BetWifi) = DetWi-ii)* (2.5) 

where the * conjugation above is understood in the lattice and 
color indices, that is U* ^ ~ ^(n+u) -v- At /i ^ the deter- 
minant is complex (while, due to the symmetries of the Yang- 
Mills integration the full partition function remains real). 



Numerical simulations are based on defining an efficient 
importance sampling of the configurations. Since the inte- 
grand (for simplicity we shall still call it 'Boltzmann factor'): 



B 



(2.6) 



is not a real, positive definite number it does not define a prob- 
ability measure for the Yang-Mills integration. There have 
been a number of methods devised to cope with this problem, 
which all involve simulating a different ensemble and correct- 
ing the results either by continuing in \x or by redefining the 
observables. 

Continuation methods use the Taylor expansion ifToll . ifTlll 
or more sophisticate expansions 11211 to enter the region of 
real, non-zero ^ by fitting the coefficients from ^ = sim- 
ulations ifioll or from simulations at imaginary [i ifllll lfl2ll . 
They rely on correctly identifying the analytic properties of 
the partition function and the various expectation values. Due 
to the noise in determining the expansion coefficients the qual- 
ity of the continuation degrades rapidly with increasing (real) 
/I. Since the simulations are done with dynamical quarks the 
statistics is limited. 

The so called 'reweighting method' proceeds by choosing a 
positive definite measure Pq obtained by splitting the original 
'Boltzmann factor' according to 



P = PqWo . 



(2.7) 



Po is used to produce an ensemble of configurations C° = 
{[/}° (where n indexes the configurations) to be reweighted 
by the complex numbers wo,n = Bn/Bo,n associated with the 
configurations C° in calculating expectation values: 



(O) 



{woO)o 
{wa)o 



(2.8) 



with O some observable and (. . . )o denoting averages over 
the ensemble C°. Notice that wq is both complex and non- 
local since it comes from the fermionic determinant. The 
(. . . )o averages contain therefore alternating contributions 
with large cancellations (the 'sign problem'). Moreover, the 
reweighting can correct an underestimated contribution in the 
C*^ ensemble, but fails if the underestimation is too drastic 
(the 'overlap problem'). In both cases the problems are ag- 
gravated by the non-locality of wq which makes it difficult to 
achieve high statistics. 

Calculations based on various implementations of the 
reweighting method ifisll have been performed mainly at small 
/i, where they agree reasonably well with other methods (an- 
alytic expansion ifTHl . lfT2ll . |0)- At large /i, however, there 
are only few numerical results yet, mainly based on only one 
method ifTsll and corroboration by different methods is miss- 
ing. 

At large /i the behaviour of QCD quantities may however 
be dominated by certain factors in the fermionic determinant 
which lead to a simpler model that is actually easier to sim- 
ulate. In its lowest order this model is considered to define 
what can be called 'quenched, non-zero density QCD' 101. 
The model is based on an analytic expansion of QCD (the 
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hopping parameter expansion) and involves the Polyakov loop 
variables of the theory, which in many setups are thought to 
catch important effects of the fermionic matter ill. This, 
and its suitability for numerical simulations makes this model 
interesting for study. Moreover it may give us hints for im- 
proving the algorithms for the full QCD at non-zero density. 

In the next subsections we shall recall the hopping parame- 
ter expansion and describe the model. 



Hopping parameter expansion of the fermionic 
determinant 



The large mass (hopping parameter) expansion of QCD 
arises from an expansion of the logarithm of the fermionic 
determinant exhibiting only closed loops: 



(2.9) 



DetW = exp(Tr InW^) 

i=l{C,}s=l 

oc 

= n n n DctD,c(i - {^f^gicc) . 

l=i {Ci} f 

Here Ci are distinguishable, non-exactly-self-repeating closed 
paths of length I and s is the number of times a loop Cci covers 
Ci . With A denoting the links along C; we have 



Cc. = 



9i 



(2.10) 

WeCi J 
f^^^±N^f.fy jf^^ ^ 'Polyakov r-path' , 
1 otherwise . (2.11) 



The index _D, C in (|2.9l l means that the traces (the deter- 
minants) are understood both over Dirac and color indices. 
A 'Polyakov r-path' closes over the lattice in the ±4 direc- 
tion with winding number r and periodic(antiperiodic) b.c. 
(e = +1(— 1)). We assume periodic b.c. in the 'spatial' di- 
rections. Notice that, since the determinant is a polynomial 
in K this expansion terminates at the order dN^Ncnf with 
d ^ 2,4 the dimension, Ni the lattice volume, Nc the num- 
ber of colors and Uf the number of flavors. For details see 

im. 



C. The massive, dense limit of the fermionic determinant 



The double limit id 



0, /i — > oo, Ke^ = C : fixed 



(2.12) 



produces a static, dense, charged background on the lattice, 
and has been therefore proposed and studied as a non-zero 
density quenched approximation |4, 5]. Note that the pure 
Yang-Mills limit corresponds to ^ = 0, which for fixed 
nonzero k requires /i — cxo. 



In the limit ( 12.12b the fermionic determinant simplifies con- 
siderably, e.g., for 1 flavor we have: 



Z™(C,{[/})-cxp 



{S} s=l 



{ecy 



Tr (Vs) 



n Det (I - e CVsf , C = (2 C)"^^ , (2.13) 



{S} 



where Vs denotes the Polyakov loop 



t=0 



(2.14) 



and from now on traces and determinants are understood only 
over the color indices. For later reference we also define the 
shortening: 



P= -^TrP. 



p* 



Nr. 



(2.15) 



(notice the different normalization to ( 12.14b above). In the 
limit ( 12.121 ) /i diverges and the parameter of the model is ( 
( 12. 12b or the related C ( 12.131 ) which is directly connected to 
the average charge density on a non-zero temperature lattice: 



no 



zfl)~2C(^Tr7'5). (2.16) 



One can study the behavior of various quantities, such as glu- 
onic correlation functions and correlation functions involving 
light quarks on such a static background, much like in the 
quenched approximation at /i = 0. However, effects expected 
to be due to the mobility of charges, in particular the possi- 
bility of new phases in dependence on the chemical potential 
cannot be studied here. 

Since this limit is obtained in an analytic expansion, we 
can systematically consider higher order corrections. In the 
following we shall study the model which is obtained at the 
next order. 



D. Large /i limit in order as a model for high density QCD 

The fermionic determinant to this order is given by: 

Zf {n, ^^, {[/}) = cxp (-2 ^ X 

r,q,i,t,t' 



X Tr 



x,r,q,i.t,t' 



(2.17) 
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in. ANALYTIC COMPUTATIONS 




K"-0(lf ) 

Determinant contributions 
for 



Contributions to order K" 



ic->0,|l->oo, 5 = Ke>'=fi,\ed C = (20"^ 



Temporal gauge 

FIG. 1 : Periodic lattice, loops, temporal gauge. In the maximal tem- 
poral gauge also the links of the basis line are fixed to I up to the 
rightmost one. 



The loops contributing to the determinant are shown in Fig. 
[T] In the following we shall use antiperiodic b.c. (e = —1) to 
ensure reflection positivity. 

For easy bookkeeping we use the temporal gauge 



UnA = 1, except for U,s„^^ 



Vs : free, (2.18) 



then 



' x,i,t,t' 



{VsY-'^u^sMVs+^^yu*.-,t'u (2-19) 

r>q>0,i^ ±1,±2,±3, 
l<t<t' <Nr (t < t' for g = 0) . 

See ^. Notice that for SU{3) we have: 



Det(I + CP) = 1 + CTtV + C^TTr* 



= 1 + 3C P + 3C^ P* + C'\ (2.20) 

Our model is thus defined by using z]^' for Z in 
Eqs. (l2.3l2.1l i rewritten for general number of flavors Uf. 

Since Zp is factorizable it is easily calculable. It is sugges- 
tive to use a splitting Eq. (12.7b preserving the factorization 
property which would allow to design a local algorithm for 
producing the C" ensemble. 

Preliminary results have been reported in 1^, ifisll . Here 
we report an extensive analysis of the phase structure of this 
model at large fi. 



A. Strong coupling/hopping parameter expansion 

As a first orientation about the behavior of the model we 
consider the strong coupling and hopping parameter expan- 
sion, which will also serve as a check of the Monte Carlo re- 
sults. For simplicity we limit ourselves to one flavor here. The 
expansion proceeds in powers of the parameters /3 and k; we 
are mainly interested in the results for the expectation values 
(Ps) of the Polyakov loop and its adjoint (Pg). 

Some details of the computation are given in Appendix A. 
The results for (P) and (P*) to order are 



1 



|C3 



1 + 4C3 

2PK^{Nr - 1) 



C6 



1 + 

2 + 3C2 + 6C6 



(1+4C3 + C6)(3 + 2C3) 



(3.1) 



and 



{P' 



\[2] 



c 



1 + 4C73 + C6 



2l3K'^{Nr-l) (1 + C3)4 + 7C6 



3 (1 + 4C3 + C6)(2 + 3C3) 

The leading behavior of this for small C is 



and 



{p*Y 



:C 



1 + -/3k2(7v^_1) )) 



(3.2) 



(3.3) 



(3.4) 



In Figs. |2]and[3]we compare the results for P and P* of the 
Monte Carlo simulations on 4^ and 6* lattices, for k = 0.12, 
one flavor and different values of /?, with pl^l and P*!^]. The 
agreement is good for the 4^ lattice and (3 — 3, while for f3 = 
5 there are already significant deviations. But the agreement 
between Monte Carlo and strong coupling results is sufficient 
to validate the simulations. On the other hand, on the 6^ 
lattice there is a remarkable difference between (3 — 5.5 and 
5.6; while in the former case the agreement with the strong 
coupling expansion remains good up to /i w 0.95 at least for 
(P), in the latter case the simulation results start deviating 
from strong coupling at much lower values of /i. This can be 
seen as an indication of a phase transition in this region. 



B. Mean field calculations 

Mean field calculations were quite popular in the early 
years of lattice gauge theory. They generally gave reasonably 
good indications of the phase structure of various models, but 



FIG. 2: Comparison with strong coupling at (3 — 3 (upper plot) and 
/3 = 5 (lower plot), 4"* lattice. Full symbols denote ReP, empty 
symbols ReP*, the lines show the corresponding strong coupling 
results. 



with the development of high speed computers and the con^e- 
sponding improvement of Monte Carlo calculations they fell 
more or less into oblivion. The reason we are reviving them 
here is to get some qualitative insight into the phase struc- 
ture of our model to which the Monte Carlo simulation can 
be compared. But it should be kept in mind that the method 
suffers from a certain amount of non-uniqueness and one has 
to apply it with some common sense. Since the mean field 
approximation of our model shows some peculiarities and has 
not been discussed anywhere in the literature, we found it nec- 
essary to derive it from the beginning. We summarize here the 
results and give details in the appendix. 

The experience with mean field theory showed that its qual- 
ity is poor without gauge fixing, but with temporal gauge fix- 
ing in pure Yang-Mills theory at zero temperature one gets 
reasonable results. Since we are dealing here with finite tem- 
perature, temporal gauge fixing is not possible. One possibil- 
ity would be the 'maximal temporal gauge' which requires to 
fix all temporal links to the identity except in one layer, but 
applying the mean field approximation would lead to a mean 
field that is not constant under time translations; this would 
not only be cumbersome, but probably also a poor approxi- 
mation since it is violating a basic symmetry of the problem. 
We take instead the next simplest choice: we fix the temporal 
gauge field to be constant ('constant temporal' or 'Polyakov 



FIG. 3: Comparison with strong coupling, P — 5.5 (upper plot) and 
(3 — 5.6 (lower plot), 6"^ lattice. Symbols as in Fig|2] 

gauge'). While the maximal temporal gauge does not lead to a 
nontrivial Faddeev-Popov determinant, going from that to the 
constant temporal gauge involves a nontrivial Jacobian (see 
appendix). 

A problem that was noted already in the eighties concerns 
the temperature dependence of the 'deconfining' phase tran- 
sition. This is not represented appropriately by the lead- 
ing mean field approximation if one uses an isotropic lat- 
tice and varies T be varying N^- We therefore fix (some- 
what arbitrarily) /3 and Nt and introduce the temperature 
through anisotropy between spatial and temporal parameters, 
see Eqs.( |2.2b .( |2.3l l. There we introduced two anisotropy pa- 
rameters 7g and jf; in principle they should both be deter- 
mined as a function of the single parameter jphys by requiring 
space-time symmetry at C = and T = 0. To leading order, 
however, we may set 7g = 7f = Iphys = 7; this is what 
was done in the computations in the appendix, since at this 
stage we cannot determine 'jphys and the mean field compu- 
tations are only meant to give a tentative picture of the phase 
structure. 

The temperature is then related to 7 by 

aT=^, (3.5) 

where the lattice spacing a is in principle determined by /?. 
(Notice that there is now a nonzero minimal temperature.) 
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Quenched QCD 
T, ~ 250 MeV 



FIG. 4: Mean field phase diagram (abscissa ji, ordinate 7 = Nt a T). 



The mean field approximation is expressed in terms of two 
different mean fields u and v for the spatial and temporal 
gauge field links, respectively. In Fig|4]we give an illustra- 
tive example, taken with /3 = 4 and N^- = 6. It shows a large 
'confinement' region for small T and ji corresponding to the 
trivial fixed point mentioned above with both mean fields u 
and V vanishing. For larger T or fi one crosses into a decon- 
fined regime with both mean fields u,v > 0. In the lower 
right corner there appears in addition an intermediate phase 
with u = 0, V > 0. The field v is close to its maximal value 1 
wherever it is nonzero, whereas u has smaller, varying values, 
depending on the region. 

Of course the fact that the mean fields u and v are ex- 
actly zero in some regions is an artifact of the mean field 
approximation; according to earlier experience already the 
next approximation in the saddle point expansion would elim- 
inate this feature. But qualitatively the mean field results in- 
dicate three phases in which different amounts of disorder 
are present: in the confined phase all the gauge fields are 
very much disordered, in the intermediate phase the Polyakov 
loops become ordered, while the spatial gauge fields remain 
disordered; finally there is the deconfined phase in which all 
the gauge fields show a high degree of order, but the Polyakov 
loops represented by v more so than the spatial gauge fields 
represented by u. In the mean field picture we present here, 
increasing 1^1 at fixed temperature, one first goes from the con- 
fined to the intermediate phase and then from there to the de- 
confined phase. This may be an artifact of the approximation 
and in reality the boundary between the intermediate and de- 
confined phases may go upward. In any case, the simulations 
to be shown in the next section suggest that by making the 
chemical potential very large at fixed temperature we end up 
in the 'half-ordered' phase. 



IV. SIMULATIONS AND RESULTS 
A. Phase diagram 

As stated in the introduction, the model we are studying 
arises from the double limit k and /i of QCD, keep- 




(infinity) 



H ~log(m) 



K~ 1 / Mass 

FIG. 5: Tentative phase diagram in T and /j, for various k. 



ing ( = K exp(/x) fixed. It can be seen either as a laboratory to 
study QCD at large mass density near the quenched limit with 
a non-zero baryon density or as a model interesting by itself 
at any value of /i and k^, describing a dense system of heavy 
baryons. 

The model still has a the sign problem that is getting more 
serious with increasing fi. But for not too large values of /i 
and not too large lattices a local algorithm with a reweighting 
still converges in reasonable computer time, as will be shown 
explicitly below. Thus we are able to carry out simulations 
across large /i "transitions" at T significantly below the de- 
confining temperature Tc at /i = 0. 

The tentative phase diagrams T vs.fi are shown in Fig. |5] 
Here we show three planes: One corresponds to "quenched" 
QCD with a finite density of infinitely heavy quarks at k = 0. 
This case has been studied for small Nr in fl Etl- At zero 
density we should find the first order phase transition of pure 
SU(3) Yang-Mills theory at « 250 MeV. 

The plane in front is the region of k near the critical value 
corresponding to masses that are small in lattice units. Here it 
has been found that there is only a crossover between confined 
and deconfined phases for all values of ^ < fXc, fJ'c ~ 400 
MeV. For ij,> Hc one expects a sharp transition, curving down 
towards T = with increasing /i [2]. It has been conjectured 
that at small T above some value of a new phase exists, dif- 
ferent from the deconfined (quark-gluon plasma) phase; this 
phase might be describable as a color superconductor and if 
the number of flavors \s N^ — i "color flavor locking" (CFL) 
is expected (O. 

Our model corresponds to a plane in between, i.e. small but 
positive K, to be chosen below; as described in Section II, it is 
based on an expansion of the hopping parameter up to order 
K^. Since k, is essentially proportional to \/M, our model 
contains some unquenched dymanics due to the fact that we 
are near but not in the quenched limit k = 0. We expect the 
phase diagram to be similar to the one for small mass just 
described. To check this is one of the purposes of this study. 

We are studying here for k = 0.12, mostly the region of 
high ji, see Fig. |6] In this region the phase diagram in tem- 
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Fixed mass plane 



this model 

(at large mass) 




y=x y=x+i 



y=x+i+j y=x 



U 



FIG. 6: Fixed mass plane phase diagram; dotted arrows indicate se- 
quences of runs. 



perature and chemical potential is expected to have a line of 
deconfinement transitions running into a triple point at some 
nonzero n and T. As mentioned above, at this point two 
further phase transition lines branch off, separating the new 
"color superconducting" or color-flavor locked phase from the 
quark-gluon plasma as well as the confined hadronic phase. It 
has been a long standing challenge for lattice QCD to explore 
this region. 



order K K 



Contributions to the quark propagator to order K ^ 



y = x 



y = x + i 



y = x 



order K K K 

Contributions to the di-quark propagator to order 
FIG. 7: Paths contributing to quark and diquark "propagators" 



B. Observables 

We measure several observables under the variation of fi 
and T, to check the properties of the different phases for small 
T and large /i. In the following we specialize to Nc = 3. The 
observables are: the Polyakov loop, 

(^) = (^Et^^.-) = (i^E^^^)' (4.1) 



and its susceptibility 



(4.2) 



where the contribution of each flavor is: 



7V3 



rib 



T3 sm 



n = ho + ni , 



no 



with the corresponding susceptibility 



XuB = ("s) - {nsf , 



(4.5) 



the (dimensionless) baryon number density hb. 



/ 



(4.3) 



the spatial and temporal plaque ttes ^Tr Paa, |Tr Po-t ^nd 

the topological susceptibility xtop = (Qtop) /i^^^r)- The 
topological charge was measured using an improved field the- 
oretical formula based on five Wilson loops ll20ll . In order 
to check the character of the conjectured third phase we also 
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measure the diquark - diquark correlators 

X J2 ([V-- t)] [^PfC^Ptiy. t + r)]*) 



-{5t6) + i5^5\){5l5t + i5^5i) 
X E W^kUi^^ i + ^)C^W-l;J{x, t-y,t + t)C 

x,y,t ^ 

-KiU^' 2/: i + y> ^ + , (4.6) 

where W^^ is the quark propagator measured in max- 
imal temporal gauge, C the charge conjugation matrix 
{a, • • • ; i, ■ • • } the color the flavor indices, respectively, and 
we have dropped the (summed over) Dirac indices. ^ is a pa- 
rameter allowing various combinations of color-flavor "lock- 
ing" (see ifigll ). Fig. |2]shows the contributions to order k^^^* 
to quark and di-quark propagators. The corresponding sus- 
ceptibility is the integral of Cqq. 

C. Algorithm and simulations 

We use the Wilson action and Wilson fermions within a 
reweighting procedure. The updating is performed with a lo- 
cal Boltzmann factor which only leads to a redefinition of the 
"rest plaquette": 

Boi{U})= n e^«^T'-^''^«x 

Plaq 

X n cxp |2Ci?eTr k- + ^ V^'l J | . (4.7) 
The weight (global, vectorizable) is 

w{{U}) = n cxp |- 2 C i?eTr Vg + E ^It.t' 

X ^ i.t.t' 



X 4^1 ({[/}), (4.8) 



such that. 



w 



Plaq 



Averages are calculated by reweighting according to 
Eqs.dfzli, dOl l. 

We have employed the Cabibbo-Marinari heat-bath proce- 
dure mixed with over-relaxation. This updating already takes 
into account part of the /i > effects and the generated en- 
semble can thus have a better overlap with the true one than 
an updating at ^ = 0. One can also use an improved Bq, 
to be taken care of by a supplementary Metropolis check. 
Anisotropy can be straightforwardly introduced. Notice that 
extracting a factor like Bq may also improve convergence of 
full QCD simulations at /i > 0. 
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FIG. 8: Data taken in the plane f3 vs. /i for fixed k — 0.12. 



The simulations are mainly done on lattice 6^ for Uf = 
1 , 3 degenerate flavors (any mixture of flavors can be imple- 
mented). The K dependence has been analyzed in |0]. Here we 
set K = 0.12 (rather "small" bare mass Mq — 0.167) which 
drives the 1/Af ^ effects in the baryonic density to about 50%. 
The task we have set to ourselves is primarily to explore the 
phase structure of the model at large chemical potential and 
"small" temperature and we accordingly vary and (3. We 
also want to check the behavior of bulk properties around the 
prospective "transition" line. 



D. Results and discussion 

The algorithm works reasonably well over a large range of 
parameters even at small temperature. The model permits to 
vary /x, k, (3 as independent parameters and it is reasonably 
cheap to measure various correlations. The region we have 
analyzed on a 6** lattice with = 3 is shown in Fig. [8] We 
have also run simulations on larger and smaller lattices, but we 
decided to base our discussion on the 6^ data and also on one 
value K = 0.12. For 8"^ x 4 and 8"' lattices the 7i/ = 3 data are 
not good enough in the (interesting) high /i region and there- 
fore we do not introduce them in the discussion. All results 
are expressed in lattice units, and we simulate the temperature 
variation by varying [3 according to ( 12.41 ) with •'jphys = 1. To 
avoid the problem of fixing the scale we shall consider T/Tc 
with Tc of the fj, = 0, pure gauge theory. We shall comment 
on all this in the conclusions. 

In Fig. |9]we show the behavior of the baryonic density ns 
with (3 at fixed /i values. We see at the different values of fi in- 
flection points (maximal slope) in /3 indicating possible qual- 
itative changes of behavior suggesting transitions from low to 
high temperature phases. In Fig. [TO]we vary fj, at several fixed 
f3 values and see the expected rapid increase of np with /i, in- 
dicating that we do not see yet saturation effects lulll . Finally, 
in Fig. [TT] we show the "landscape" of the real part of the 
baryon density (while the imaginary part is compatible with 
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FIG. 9: Baryonic density vs. j3 at fixed ji. 
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FIG. 10: Baryonic density vs. /i at fixed /3. 



zero inside the statistical errors, as it should be). 

A clearer view of the situation is provided by looking at the 
"landscape" of the susceptibility of the baryon density, which 
is shown in Fig. [T2] A ridge is clearly visible, highlighted by 
a dashed black line. A second line (dotted) will be explained 
later. 

The main variation in the baryon density is an exponen- 



tial growth with ^. This masks to a certain extent the finer 
structure. We found it therefore advantageous to look at the 
Polyakov loops and their susceptibility. In Fig. [T3]we show 
this susceptibility at fixed /i vs (3 and in Fig. [T4]at fixed (3 vs. 
/i, and in Figs. [15] and [16] the corresponding landscape. 

The plots of the Polyakov susceptibility show quite clearly 
maxima indicating possible transitions or crossovers. In the 
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FIG. 11: Landscape of the baryonic density. The color scale (right) 
is based on logj^g (ns ) . 
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FIG. 12: Landscape of the baryon density susceptibility. The color 
scale (right) is based on logj^,) (xn^ ) • 



landscape Figs. [T5]and[T6lone of these maxima shows up as a 
well defined ridge, indicated by a dashed black line. It shows 
only a moderate slope in i^i, which explains why the maxima 
are more pronounced when we vary f3 at fixed /i than vice 
versa. The broadening of this ridge at small as well as of 
the maximum in Fig. [T3]is responsible for the loss of a sharp 
transition signal at small ji. These figures clearly show that 
the transition at fixed ji ~ 0.50 is less steep than the one at 
ji = 0.80. Presumably at /i <~ 0.6 we are dealing with a 
crossover, whereas at large /^t the signal is more compatible 
with a real phase transition. Notice that changing /3 at fixed /i, 
we cross the transition line at a more oblique angle at smaller 
/i, but the broadening of the ridge and loss of a transition sig- 
nal is a genuine effect, as can be seen from Figs. [T5]and[T6l 

A second ridge branching off from this main ridge at large 
/i, highlighted by a dotted line is suggested by looking at the 



level lines in Fig. [T5]and corresponds to the second maximum 
suggested at large /^t in Fig. [14] This may indicate the appear- 
ance of the new phase at large /i and small T/Tc discussed 
above. 

We use the results for the Polyakov loop susceptibility to 
estimate the possible position of the transition points in the /3 
vs /i plane; to go half way toward a possible physical interpre- 
tation the positions determined in this way are indicated by 
the blobs in the diagram T/Tc vs. fiphys/Tc of Fig. [17] where 
fJ-phys ~ fJ-/o,iP) = NrnT and the relation between j3 and 
T/Tc has been roughly estimated from the // = quenched 
QCD with N-r = 6 (we shall comment on this point in the 
conclusion section). In this figure the axis of the blobs indi- 
cate the search lines in the simulation. The shaded blobs cor- 
respond to the rather unambiguous 'deconfining' signal ob- 
served for fj. >~ 0.6 {f3 <^ 5.72). The 'transition' line 
suggested by this signal starts at the lower point A on the 
figure, located at /? ~ 5.55, /i ~ 0.88, i.e., with our rough 
estimation fiphys/Tc — 2.4, T/Tc — 0.45 (below which we 
could no longer obtain reliable data) and ends at the point B 
located near (3 ~ 5.72, jj, ~ 0.6, i.e., with our rough esti- 
mation ^iphys/Tc — 2.3, T/Tc — 0.65. Above this point the 
signal becomes ambiguous. But one should keep in mind that 
moving along lines of fixed ji across a broad ridge, the max- 
imum in general is shifted with respect to the ridge (in our 
case to lower j3 values), the location of a transition becomes 
somewhat blurred, in accordance with the claim that here we 
are dealing with a crossover and not a phase transition. In Fig. 
[TTlwe shaded the upper, 'broad ridge region' above B where 
the maximum at fixed /i or /3 deviates significantly from the 
location of the ridge, which can be easily understood from the 
landscape Fig. [15] Notice that since we keep k fixed /i = 
does not represent the pure Yang Mills theory therefore we did 
not try to go to this limit. The white blobs correspond to the 
more volatile, possible 'transition' branching off near point A 
at large /i, whose signal is strongly affected by fluctuations. 
We also shaded the region at high /i in the lower right hand 
corner, where we could not obtain reliable data due to the sign 
problem. 

The picture emerging from the data is thus the following: 
for ji < 0.5 — 0.6 iiJiphys/T ~ 3) there is only a broad 
crossover, while for 0.6 < ^ < 0.9 (3.6 < fiphys/T < 5.3) 
there is evidence of a sharper crossover or transition at a value 
l^c depending on /3. Moreover, for jj, ~ 0.9 there is some ev- 
idence of the presence of the second transition even though 
this evidence is much weaker than the other one because at 
larger values of /i the fermion determinant strongly oscillates 
and, indeed, the usual sign problem manifest its effects. 

To get some further insight into the nature of the different 
regimes or phases we also wanted to look at the distribution 
of the values of the Polyakov loop in the complex plane. At 
first we considered the 'histograms' corresponding to the fol- 
lowing mathematical expression; 

HA{x,y) = 
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where x is any point in the spatial lattice and 0a.s(^) is the 
function which is 1 if |t — .s| < A/2 and otherwise (the 
arguments x,ij in H should not be confounded with space- 
time points). For the figures we used 20 x 20 bins choosing 
A accordingly. These quantities have the advantage that they 
are positive, because they use the expectation values (.}o de- 
termined by the positive Boltzmann factor Bq (see Eq. 12. 7t ; 



therefore they can be interpreted as probability distributions. 
But their disadvantage is that they depend on the choice of 
Bq. It should also be noted that they describe not really the 
distribution of the Polyakov loops themselves, but rather the 
product of the Polyakov loop with the weight factor w; for this 
reason absolute values larger than 1 are possible and actually 
occur, as we will see. 
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FIG. 15: Landscape of the Polyakov loop susceptibility. The color 
scale (left) is based on logjg(xp) 
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FIG. 16: 3d view of Fig. [B] 

As an example, see Fig. [18] and Fig. [T9]that represent the 
histogram of Ha at different values of /i at /3 = 5.65 and 
different values of /3 at /i = 0.70, respectively. These figures 
show different behavior of this observable in accordance with 
the transition lines indicated in Fig. [17] In fact in Fig. [T8]one 
can discern three different regions: the first one corresponds 
to n < 0.6, where the Polyakov loops are concentrated in 
a small region around zero with only a slight preference for 
positive real parts; in the second region, for 0.6 < fi < 0.9 the 
Polyakov loops become considerably larger, favoring positive 
real parts in a significant way, while finally for /i > 0.9 the 
Polyakov loops (times weight) becomes quite large, but are 
distributed almost symmetrically around the origin. 

This picture can be corroborated by looking at Fig. [19] 
which according to Fig. [17] should only show one transition. 
One can see a change of behavior around the point (3 — 5.65 
(which also occurs in Fig. [TSj: The Polyakov loops become 
somewhat larger with a distribution more heavily favoring 



positive real parts; we interpret this as the transition from a 
confined to a deconfined phase. 

A 'distribution' independent of the choice of Bq can be de- 
fined by considering 

Ta(x,2/) = {eAARePs)eA,yiImPs)) , (4.10) 

which means adding the weights of all configurations produc- 
ing a value in a given bin — a;| < A/2 , \RePg — 
y\ < A/2. Because now the "expectation value" (.) refers to 
the complex "Boltzmann factor" B (see Eq. 12.7b . Ta is com- 
plex and does not represent a probability distribution. But for 
small A we have 

{P)^Y.^x + ty)TA{x,y), (4.11) 

where the sum runs over a lattice with lattice constant A in 
the xy-plane. Since the expectation value of P is real, RcTa 
has to be even and IttiTa odd in y. 

We give some representative figures showing the behavior 
of Ta across the putative transitions, for the same parameters 
as before. Fig. [20l shows RcTa for /3 = 5.65 for various in- 
creasing values of ^. Again we should observe the crossing of 
two of the putative transition lines. The transition signals are 
not very strong, but we can observe that for /.t < 0.7 negative 
real parts are present, which disappear for > 0.7; at /i > 0.9 
the real parts become considerably larger again, reaching val- 
ues of 0.3. Fig. l2n shows ReT at n — 0.7 for increasing val- 
ues of /3. Here the parameters are such that we should observe 
only the transition between the hadronic and plasma phases. 
The indication for this is again that the real parts touch the 
origin for /3 < 5.65, whereas for /3 > 5.65 they increase to 
positive values, but staying below 0.2. 

Both Fig. |20]and Fig. 1211 show that RbTa is to good accu- 
racy even in y, as required for the reality of {P). In Figs [22] 
and|23]we show the imaginary parts of the 'distributions' Ta. 
The qualitative signal of the transitions/crossovers is similar 
to that of ReT A- It should be noted that now IitiTa is, to 
very good precision, odd in y, again in agreement with the 
reality of (P). 

Polyakov loops and charge density (and their susceptibil- 
ities), have been the primary quantities used to uncover the 
phase structure. We also have measured plaquette averages 
(for both temporal and spatial plaquettes), the topological 
charge density (using the improved field definition) and quark 
and di-quark correlators (in maximal axial gauge). All these 
quantities also some show peculiar behavior in both n and /3 
which will be exemplified here on two chosen runs, at fixed 
(3 ~ 5.65 vs. /i and at fixed /i ~ 0.7 vs. (3: In Figs. [24] 
and [25] we present the dependence of the plaquette averages 
on /I at /3 = 5.65 and on /3 at /i = 0.7, respectively. We 
see here clearly the emergence of a physical energy density 
by the gap developing between the spatial and temporal pla- 
quettes with increasing /i and (3; this corroborates the phase 
picture derived before. In Figs. [26] and [27] we present for 
the same runs the topological susceptibility whose behavior 
again is in agreement with the previous conclusions since it 
decreases in the region where we expect deconfining to set in. 
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^lphys/Tc = 6^T/Te 



FIG. 17: Phase diagram in the P (or T /Tc) - f^phys/Tc QCD plane. The dotted straight lines correspond to constant /i, the dashed ones to 
constant /3. The blobs, shadowing and other features are explained in the text. 



Finally in Figs. |28]and|29]we present the dependence on /i and 
on P of the diquark susceptibility obtained by integrating the 
diquark-correlators Eg. (14. 6b for ^ = 0.5; here we only show 
the contribution to this susceptibility from the terms. This 



corresponds to quarks showing a (limited) amount of mobil- 
ity and as can be seen from these figures, the susceptibility 
to this order is sensitive to the chemical potential (while the 
zero-th order contribution is dominated by a contact term and 
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FIG. 18: Polyakov loop 'histogram' Hj:^{x, y) of eq. ( 14.9b vs. /i at /? = 5.65. 




FIG. 19: Polyakov loop 'histogram' H/\{x, y) of eq. ( |4.9l l vs. /3 at /i = 0.70. 
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FIG. 21: Real part of the Polyakov loop 'distribution' Ti^{x,y) of eq. l l4.10t vs. /3 at /i = 0.70 fixed. 
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FIG. 23: Imaginary part of the Polyakov loop 'distribution' Ta{x, y) of eq. | |4. lOt vs. /3 at jj, — 0.70 fixed. 
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FIG. 24: Plaquette averages vs. fi at fixed /3 — 5.65. 
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FIG. 26: Topological susceptibility average vs. /i at fixed /3 — 5.65. 
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FIG. 25: Plaquette averages vs. /3 at fixed /j, = 0.70. 



FIG. 27: Topological susceptibility average vs. (3 at fixed /i — 0.70. 



is rather flat). The strong increase with ii, compared with the 
rather flat /3 dependence may indicate new properties of the 
matter at high density. 



V. CONCLUSIONS 

To obtain analytic informations about our model we first an- 
alyzed it via the strong coupling expansion; the agreement for 
/3 < 5.5 and small /i with the numerical simulations should 
be seen as a validation of the simulation program. But our 
calculations show strong effects at slightly larger ^ , which 
akeady at /3 = 5.6 depart considerably from strong coupling 
estimates; this is an indication of a possible phase transition. 
Next we obtained a phase diagram in a mean field approxima- 
tion, showing the existence of three different phases. 

The phase structure found by the numerical simulations for 
Uf = 3 is shown in Fig. [17] The signal for the deconfining 
transition (or narrow crossover) on the line connecting A and 
B is rather good and it also appears that at small /i (above B) 
the transition is smoothed out in accordance with the expecta- 
tions from full QCD simulations ll2ll. ll22ll . A second transition 



at large i^i could only be identified tentatively. In this region, 
the diquark susceptibility grows strongly. This region needs 
further study to reach a conclusion, but it is interesting that 
the general picture shows qualitative agreement with the one 
found in the mean field approximation. 

The algorithm works reasonably well over a wide range of 
parameters and for lattices up to 6^ (S"' for Uf = 1). We obtain 
large densities for temperatures ^ ^ Tc or less and reach ratios 
^p^"^^ 5. It appears difficult, however, to go to larger lat- 
tices and larger ^ with this algorithm and one should consider 
improving it. For the time being, however, these difficulties 
precluded us from performing further tests, such as finite size 
analysis, in order to establish unequivocally the character of 
the various transitions. 

The model permits to vary /i, k, (3 and Nr as independent 
parameters. Also anisotropic lattices can be envisaged. It is 
therefore interesting to extend the study to take advantage of 
this full variability. Also extending the model to higher or- 
ders in K can be envisaged. The bookkeeping soon becomes 
unmanageable, one could however consider using statistical 
ensembles of large loops ll23ll . 

A related matter is the relation to physical quantities such 
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FIG. 28: Diquark susceptibility average vs. /i at fixed (3 = 5.65. 



gluon dynamics of the SU(3) theory in this situation. It would 
then be natural to think of it as providing a heavy, dense, 
charged background for propagation of light quarks and calcu- 
late light hadron spectra and other hadronic properties under 
such conditions. This could also help fixing a scale control- 
ling the behavior of the light matter We consider pursuing 
work on this subject. 
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FIG. 29: Diquark susceptibility average vs. /3 at fixed jj, = 0.70. 



as temperature and masses. In this study we introduced a T- 
dependence by varying (3 and tried to avoid the necessity of 
defining a scale by considering only dimensionless ratios such 
as Uphys/T. This, however, has to be taken with a grain of 
salt: indeed, varying (3 also introduces varying finite volume 
and quark 'mass' effects. It would be less ambiguous to vary 
Nr if we could reach large enough lattices. Alternatively one 
could consider using a variable anisotropy. In a first approxi- 
mation one could take 7^ = 7f = Jphys, such as in the mean 
field approximation in section III.B, but non-perturbative cor- 
rections might be large and a bona-fide calibration may be- 
come necessary |@]. All renormalization questions, however, 
are difficult when we need to consider the effects of the quarks 
as introduced in fixed order hopping parameter expansion. 

Concerning the significance of this analysis we can take two 
points of view: 

Firstly, we can consider this model for itself, as describing 
'quasi-static charges' interacting via gauge forces and having 
a non-trivial phase structure. 

Secondly, we can consider this model as an evolved 
'quenched approximation' in the presence of charged matter. 
Then this study would give us information about the modified 



1. Strong coupling expansion: some details 

We first calculate the term of order zero, which would van- 
ish trivially without the presence of the chemical potential 
term C. The fermion determinant to order k° is 



2P ^lldct{l + CVs)\ 



(A.l) 



where the determinant only refers to the color degrees of free- 
dom. In order to evaluate this explicitly we introduce the char- 
acters X(j of the irreducible representations a of SU (3). In the 
maximal temporal gauge Vg is simply given by Vg and we find 

ZP = n (1 + CXiiVs) + C\-3{Vs) + 6-3)'. (A.2) 

X 

Using the well-known facts (see for instance ll24ll25ll ') 



X3X3 = Xi + X8 , 

X3X3 = XS + X6 , 
X3X3 = X3 + Xg , 

and defining D = 1 + 4C^ + this becomes 

2C + 



(A.3) 



3C2 + 20"^ 



D 



1 



D 



-XsiVs) + ^C\e{Vs) 



(A.4) 



From this it is straightforward to obtain the expectation values 
{P3) and to order as 



{P) 



[0] 



1 



|C3 



1 + 4C3 + 



and 



' ' 1 I / 



C3 



1 + 4C3 + C6 



(A.5) 



(A.6) 
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The next nontrivial order is O(k^) in the fermion determinant 
and comes from the Polyakov loops with one excursion to a 
neighboring site. A nonzero result is obtained only by com- 
bining it with terms from the Yang-Mills action; the lowest 
nontrivial contribution is therefore 0(k^/3). Concretely we 
obtain to order 



2 [2] 



/ J ' x,i.t,t' 

x,i.t,t' 



(A.7) 



After integrating over the spatial gauge fields U only terms 
with t' = t + 1 survive; the integrals occurring are of the form 

(A.8) 

Thus we obtain before the integration over the V's 

J n dUZP = ZP + f3Cx3m^ (A.9) 

withC = 2/3C{Nt- — 1)k'^ /3- To obtain the expectation values 
of the Polyakov loops from this we have to expand the product 
in irreducible characters; we need only the terms involving the 
representations 3, 3, 1. Using Eq.( IA.4b we see that we need a 
few more decompositions of SU{3) representations, namely 



X3X6 

XaXe 

X3X8 



X8 + Xio 

X3 + Xl5 

X3 + X6+ Xl5 



(A. 10) 



Since the expectation values are normalized by the partition 
function, as usual only connected contributions occur; thus 
the results for {P) and (P*) to order are 



1 



|C3 



1- 



1 + 4C3 + C6 



(1 + 4C3 + C6)(3 + 2C3) 



(A. 11) 



and 



C3 



1- 



1 + 4C3 + C6 

2f3K'^{Nr - 1) (1 + C^y + 7C6 

3 (1 + 4C3 + C6)(2 + 3C3) 

We note the leading behavior for small C: 



(A. 12) 



2. Mean Field: some details 

We first compute the Faddeev-Popov determinant J{v) for 
the Polyakov gauge, which can be computed as the Jaco- 
bian for the transformation from the maximal temporal to the 
Polyakov gauge. 

The reduced Haar measure for the conjugacy classes [U] of 
SU{N) is given by iH] 



sm 



dcpi . . . d(j)N-i , (A. 15) 



where A/^ is a normalization constant; this would be the appro- 
priate measure for the temporal gauge field in the unfixed links 
of the maximal temporal gauge. We are instead spreading the 
field uniformly over Nr links such that we want to integrate 
over V e SU{N) with 1/^^ = U, so we want to write 



d[U] = J{V)d[V] , 



(A. 16) 



where J{V) is now the 'quotient' of the Haar measures for 
V^^ and U, i.e. 



J{v) = -[['■ 



sm 



i<j 



2 ^ jV^(0i-0,) ^ 



2 / 



(A.17) 



So we have to integrate the homogeneous temporal gauge 
fields with the measure 



d[V] = Yl sin^ 



i<3 



Nr{4>^-(|)J) 



N-l 

fc=i 



(A. 18) 



The range of integration is the interval [— tt, tt) for each 0^; 
this means of course that covers the group SU(^]V) JVr 
times; this 'over-counting' is necessary, since otherwise the 
integration of functions of V would involve some completely 
arbitrary choice of the 'A^^-th' root. 

We now proceed in the standard fashion to produce the 
mean field theory as a saddle point approximation (see for in- 
stance ||27[ I28I1 ) for the partition function: first the integrals 
over the group SU{N) are replaced by integrals over the em- 
bedding matrix space Mpf^pf{C) by inserting the identities 



1 = J duS{U-u) 



J dM j exp [ii?eTrAf1'([/-M)];A.19) 



and 



p[2] ^ ( 1 + -pK^iNr - 1) 



(A. 13) 



(A. 14) 



for each spatial link and similarly for V, introducing the ma- 
trix valued fields v and K for each temporal link. The group 
integrals for the different Unks are then decoupled and reduce 
to the one-link integrals 



df/exp(i?eTrA/t{7) 



(A.20) 
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and 

J dVJ{V)cxp{ReTrK^V) . (A.21) 

Carrying out the integrals over the gauge field, using these 
definitions, the partition function reduces to an integral over 
the matrix valued fields u,v,M,K with an effective action 
S{u, V, M, K). This integral is suitable for a saddle point ap- 
proximation. By symmetry there must be a translation invari- 
ant extremal of S. For the matrix valued fields we further- 
more make the ansatz that they are multiples of the identity; 
by slightly abusive notation 

u = wl, ?j = (vi + 1^2)1, 
M = (mi + 17712) I, K = (fci + 7fc2)I ■ (A.22) 

We anticipated here already that u will be real. Using this 
ansatz and introducing a single asymmetry parameter 7 = 
IG = If, as discussed in Section |III] the action per site s 
becomes 

- s =2,-u^ + ?,(3-iu^{vl+vl) 

+ 3 lnC(«77ii) + ln77(iA:i, 7fc2) 

— 37(fciWi + k2V2) ~ 9imiu (A. 23) 

where the functions ( and r] are defined for arbitrary complex 
arguments a, bi, 62 as 



CO 



and 



a) = J d[U]exp{aReTTU) (A.24) 



»i,62) = J d[V]J{V)exp{biReTTV + b2 ImTiV) 

(A.25) 

the functions C, and 77 in more 



For the group SU{3) we write 
explicit form: 

d(j)i / d02Pl (</>!, 02) 

-TT J —TT 



X 



cxp [a(cos 01 + cos 02 + cos(0i + ^2)] (A. 26) 



and 



d<f)l / d(j)2PN^{4>l-,4>2) 

-TT J —Tl 

X exp [61 (cos 01 + cos 02 + cos(0i 
exp [62 (sin 0i + sin 02 — sin(0i H 



X 



'02)] 
02)] , 

(A.27) 



with 



Pfc(01, 02) 



. 2 /^M^l-02)^ • 2/^fc(01 +202) 

sm sm ' 



V 2 

2 //c(02 +20i) 



(A.28) 



V 2 

When searching for a saddle point we have to allow all pa- 
rameters to be complex. The saddle point equations, requir- 
ing stationarity of s with respect to u,vi^V2,a — mii, hi = 
iki, 62 = *^"2 are 



4/3 



bi - 



62 - 



V2 



37 3' 

ACK^{Nr - l)u{vi + ^7^2)^" , 
2 

Z — 7i fl 

7 

2C(7;i + iv2)^^'\l + 3Nr{Nr - 1)k^u'^) , 

2 — 7t 7^2 

7 

2iC(7)i + i7;2)^^"^(l + 3Nr{Nr - 1)^27*2) , 

i^ln„((.i.6.), 

1 5 

-— In 77(61,62) . (A.29) 

3 ao2 



The system of equations is of the form of a fixed point condi- 
tion and is solved by iteration. There is always a trivial fixed 
point 



77 = 7)1 = 7;2 = a = ai = ^1 = 



(A. 30) 



In general if there is more than one fixed point (which may be 
reached by choosing different starting points for the iteration). 
It turns out that all the fixed points satisfy a = 77771 real, &i = 
7fci purely imaginary, 7J2 = and 7t, vi real; note that 712 = 
is consistent with these equations because of the symmetry 



77(61, 62) = 77(61, -62), (A.31) 



which follows from the unimodularity {d[U] = (i[C/^]). 

With our sign convention one always has to choose the fixed 
point leading to the highest value of the free energy density 
f = s for the parameters chosen. This leads to discontinuities 
in the first derivative, typical for first order phase transitions, 
and finally to the phase diagram shown in Fig|4] 
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